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ABSTRACT 

We present a detailed description of a phenomenological H2 formation model and local star formation pre- 
scription based on the density of molecular (rather than total) gas. Such approach allows us to avoid the 
arbitrary density and temperature thresholds typically used in star formation recipes. We present results of the 
model based on realistic cosmological simulations of high-z galaxy formation for a grid of numerical models 
with varied dust-to-gas ratios and interstellar far UV (FUV) fluxes. Our results show that both the atomic-to- 
molecular transition on small, ~ 10 pc scales and the Kennicutt-Schmidt (KS) relation on ~ kpc scales are 
sensititive to the dust-to-gas ratio and the FUV flux. The atomic-to-molecular transition as a function of gas 
density or column density has a large scatter but is rather sharp and shifts to higher densities with decreasing 
dust-to-gas ratio and/or increasing FUV flux. Consequently, star formation is concentrated to higher gas surface 
density regions, resulting in steeper slope and lower amplitude of the KS relation at a given Eh, in less dusty 
and/or higher FUV flux environments. These trends should have a particularly strong effect on the evolution of 
low-mass, low surface brightness galaxies which typically have low dust content and anemic star formation, but 
are also likely to be important for evolution of the Milky Way-sized systems. We parameterize the dependen- 
cies observed in our simulations in convenient fitting formulae, which can be used to model the dependence of 
the KS relation on the dust-to-gas ratio and FUV flux in semi-analytic models and in cosmological simulations 
that do not include radiative transfer and H2 formation. 

Subject headings: cosmology: theory - galaxies; evolution - galaxies: formation - stars;formation - methods: 
numerical 



1. INTRODUCTION 

Conversion of gas into stars is one of the major sources of 
uncertainty in modeling formation of galaxies. This uncer- 
tainty reflects our incomplete understanding of the process of 
star formation both locally and on global scales. Traditionally, 
star formation is included in cosmological simulations and 
simulations of isolated galaxies by using simple phenomeno- 
logical prescriptions that relate local rate of star formation to 
the local density of gas, with some additional criteria such as 
temperature and density thresholds for the gas to be eligible 
for star formation. The parameters of these prescriptions are 
chosen so that the empirical power law relation between the 
surface density of star formation, Ssfr, and surface density of 
(hydrogen) gas av eraged o n kpc scales, '^sfr with 
n * 1 - L4, (Schmidt 1959; Kennicutt 1998; Bigiel et al.|2008]l 
observed in z ^ galaxies is reproduced (see, e.g., Schaye & 



|Dalla Vecchia|2008] for a recent overview). 

However, both theoretical considerations and observational 
evidence indicate that such approach may miss some impor- 
tant environmental trends. For example, relation between 
the local star formation recipe and the large-scale Kennicutt- 
Schmidt (KS) relation is not trivial and depends on the den- 
sity and therma l structure o f the interst ellar medium (ISM, 
Kravtsov|200"3i j Tassis 20 ()7i |Wada & No rman 2007; Robert- 
son & Kravtsov 2008(|Schaye & Dalla Vecchia|2008( |Saitoh| 



et al. 2008). This is because for a given large-scale gas surface 
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density the fraction of dense, star forming gas is determined 
by the gas density distribution functi on, which, in turn, de- 
pends on the thermal state of the ISM ( |Wada & Norman|20OTl 
Robertson & Kravtsov 2008|l. For the same reason, the global 



rate of star formation may be controlled by the rate with which 
dense gas is formed by the ISM, rather then by the assumed 
local efficiency of the gas (Saitoh et al. 2008). This implies 
that star formation parameters tuned to reproduce the empir- 
ical KS relation in o ne situation (e.g., in controlled simula- 
tions of isol ated disks Sp ringel & Hernquist|2003| Schaye & 
[Dalla Vecch ia 2008) may not reproduce this relation in galax- 
ies with significantly different ISM density distributions. 

In addition, there is a growing observational evidence that 
the KS relation is more complex than previously though t 
( Heyer et al.|[2004} [Boissier et aL][2003t [Bigiel et al.|[2008l ). 
For example, mstead of a well-defined surface density thresh- 
old at low £h below which Ssfr drops to zero ( [Martin & Ken-| 
|nicutt[[200T] l, observations indicate continuous rel ation be- 
tween star formation rate and gas surface densities ( Boissier 
2007|) down to sm all £h, albeit with a steeper slope 
Bigiel eraL]|2008l l 



et al. 



dwarf 



Likewise, studies of individual 

galaxies, which typically have low gas surface densi- 
ties (Eh ^ 10-20 Mq pc"^) throughout their disks, show that 
the KS relation in such galaxies is generally characterized by 
a considerably steeper slope, n =a 2-4, than the canonical 



value of 1.4 ( [Heyer et al.][2004) [Bigiel et al.[[2()08) [Verley 
et al. 2010[ l. Moreover, recent detailed study of the global 



star formation relation by Bigiel et al. ( 2008 1 shows that a 
single power law is in general a poor description of the KS 
relation over the entire range of surface densities. Instead, 
the slope of the EjpR - Eh relation may vary from the steep 
values of 2 - 4 at Eh < lOMopc"^ to linear n a; 1 at 
Eh ~ 10 - 100 Mq pc"^ and then possibly steepening again to 
« a; 1.5 - 2 at Eh > lOOMopc^^. 
Finally, the growing evidence indicates that in high-redshift 
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galaxies (z > 3) the KS relation is significantly steeper and has 
an order of magnitude lower amplitu de at Esfr ^100 Mq pc'^ 
(|Wolfe & Chen' 2006||Rafelski|2009| see also Fig. 3 in |Gnedin 
[& Kravtsov 2010^"^ 

This complex behavior of the star formation rate density 
with the density of the neutral gas (H I+H2) can be understood 



if star formation occurs only in the molecular gas (p^obert- 
son & Kravtsov"2008', 'Gnedin et al."2 009t |Krumholz et al. 



2009b ; Pelupessy & Papadopoulos 2009| Gnedin & Kravtsov 
^10) 1. Indeed, detailed observations of nearby galaxies show 



that star form ation corre lates most strongly with the molec- 
ular gas (e.g., [Wong & B litz 2002a; Bigiel et al. 2008]l, es- 



pecially with the densest gas traced by HCN emission ( Gao 



& Solomon 2004 Wu et al. 2005 1, while it only correlates 
weakly, iTat all, with the density of atomic gas (Wong &" 



Blitz||2002a| [Kennicutt e t al. 2007, Bigiel et al. 2008)1. We 



can thus expect that the relationship between the star forma- 
tion rate density and gas density Zh - + (the KS re- 
lation) varies depending on the molecular fraction of the gas 

Several factors may control the molecular fraction in the 
gas on different spatial scales. On small scales of indi- 
vidual molecular complexes it is primarily the cosmic dust 
abundance and the interstellar FUV radiation that control 



the atomic-to-molecul ar tr ansition (e.g ., Elmeg reen 1993 
see 



Krumholz et al. ||200g 

gogical review) 



Stabl er & P alla 2005TTot pHa- 
kpc) scales the fraction of 



On larger ( 

dense, molecular gas in a patch of gas of a given Xh is ex 
pected to depe nd on the den sity distribution of gas in that 
patch (e.g., ,Elmegreen||200 2|. The density dist ribution itself 



depends on the rmodynamics of gas (see, e.g., Robertson & 
Kravtsov|2008| l and metallicity, as more metal rich gas may 



be more efficient in building regions of higher densities via 
radiative shocks arising in the highly turbulent medium of 
gaseous disks. The density PDF should also reflect the global 
dynamics of gas in galactic disks in general. For example, 
spiral density wave will compress the gas facilitating its cool- 
ing and conversion of atomic gas into molecular form. Like- 
wise, large-scale instabilities seed the turbulence in the disk 
that can shape the global density PDF ( Wada & Norman 200 Ij 
[Elmegreen 2002; Kravtsov 2003, K rumholz & McKee 2005]. 

Although observational studies of environmental depen- 
dence of the KS relation on gas metallicity, interstellar FUV 
radiation, and other properties of galaxies are in their early 



sta ges (e g., [Bigiel e t al. 2008; Krumholz et al. 2009a ; |RafeF 
|ski||2009| , it is clear that such strong dependences can have 
important implications for our understanding of ga laxy evo- 
lution (see discussion in Gne din & Kravtso"v[|201()| l. For ex- 
ample, given that observations indicate that star formation in 
low-metallicity, high-UV flux environments of high-redshift 
galaxies is concentrated to significantly higher gas surface 
densities (Wolfe & Chen 2006, Rafelski 2009), stars in these 
galaxies should be confined to the high surface density re- 
gions and should therefore be more resistant against dynam- 
ical heating in mergers. At the same time, the longer gas 
consumption time scales in lower density regions of high-z 
gaseous disks along with high accretion rate would keep them 



gas rich and more resilient to mergers as well (e.g.. Robertson 
let al.|2004 2006J||Springel & Hernqmst(200 5 ). This can help 
to resolve one of the major puzzles of hierarchical galaxy for- 
mation: prevalence of thin disks at low redshifts in the face of 
high merger rates at high redshifts. 

It is thus important to explore potential effects and impli- 
cations of the enviromental dependence of the KS relation for 



the evolution of galaxies. However, to capture the key physics 
responsible for this dependence in cosmological simulations 
of galaxy formation is challenging, because this requires high 
spatial resolution to model dynamics of interstellar medium 
in the hierarchically forming galaxies, 3D radiative transfer 
to model local UV radiation flux, and formation of molecular 
hydrogen. The latter is mediated by dust grains which cat- 
alyze H2 formation and provide the initial key shielding from 
interstellar FUV radiation. This shielding allows build-up of 
molecular fraction sufficient for H2 self-shielding, which in 
turn shapes the sharp transition of atomic to molecular gas. 

Although fully self-consistent modeling of dust chemistry 
and H2 formation is still far beyond reach, phenomenolog- 
ical model capturing the essential metallicity and UV flux 
dependence of molecular fraction can be used to model H2 
in self-consistent, high-resolution cosmological simulations 
( [Gnedin et a l. 2009; Gnedin & Kravtsov 2010). In this study 
we present a detailed description of such H2 formation model 
and local star formation prescription based on the density of 
molecular (rather than total) gas. We present results for a 
grid of numerical models with varied dust-to-gas ratios and 
interstellar FUV radiation fluxes and explore the dependence 
of atomic-to-molecular transition on small, molecular cloud 
scales, on these variables and the effect this dependence has 
on the Kennicutt-Schmidt relation on large ~ kpc scales. We 
parameterize the dependencies observed in our simulations 
in convenient fitting formulae, which can be used to model 
the metallicity and UV flux dependence of the KS relation in 
semi-analytic models and in cosmological simulations that do 
not include radiative transfer and H2 formation. 

2. SIMULATIONS 

For our te sts we use the simul ation of galaxy formation 
described in [Gnedin et al.| ( 2009| l. The simulation was run 
with Adaptive Refinement Tree (ART) c ode (Kravtsov 1999j 
[Kravtsov et~aL] |2002; Rudd et al. 2008) and follows a La- 
grangian region corresponding to five virial radii of a sys- 
tem, which evolves into a typical halo of an L» galaxy (M x 
10'^ Mq) at z = 0. The mass resolution in the high-resolution 
Lagrangian region is 1.3 x 10^ Mq in dark matter and mass 
resolution in baryon that varies from ~ 10^ Mq to ~ 10^ Mq 
depending on the cell size and density. The simulation reaches 
peak spatial resolution of 260 comoving pc (65 pc in physical 
units at z = 3). The Lagrangian region is embedded into a 
cubic volume of 6/1"' comoving Mpc on a side to model the 
tidal forces from the surrounding structures properly, but this 
outer region is resolved only coarsely with a uniform 64^ grid. 

The cosmological simulation follows collapse of dark mat- 
ter and gas self-consistently. The heating and cooling of gas is 
followed as well, so that gas can dissipate the energy it gains 
during collapse and sink to the center of its parent halo. Our 
simulations include 3D radiative transfer (RT) of UV radia- 
tion from individual stellar particles formed during th e course 
of the simulation using the OTVET approximation (Gnedin 



& Abel| 2001 1. Inclusion of the RT is important because the 
local U V flux can set ionization and heating balance of gas 
and influence the abundance of molecular hydrogen, as we 
descibe below and in the Appendix. Unlike the IGM after 
reionization, which can be assumed optically thin to ionizing 
radiation, the dense ISM gas of simulated galaxies may well 
be opaque to ionizing photons of all but the nearest stars. 

The simulations incorporate non-equilibrium chemical net- 
work of hydrogen and helium and non-equilibrium cooling 
and heating rates, which make use of the local abundance of 
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atomic, molecular, and ionic species and UV intensity. This 
network includes formation of molecular hydrogen both in 
the primordial phase and on dust grains. The abundances of 
the relevant atomic and molecular species are therefore fol- 
lowed self-consistently during the course of the simulation. 
The heating and cooling terms in the equation for the inter- 
nal energy include all of the terms normally included in the 
simulations of first stars and in the ISM models, including 
cooling on metals. We describe all included reactions and 
heating/cooling processes in Appendix. 

The model also accounts both for self-shielding of H2 from 
the dissociating FUV radiation and the shielding provided by 
the interstellar dust using phenomenological prescriptions for 
shielding factors. The details of the model are presented in 
the Appendix. Our model is calibrated against the observed 
column density dependence of atomic and molecular gas frac- 
tions in the Milky Way, LMC, and SMC (see Appendix). In 
particular, the model reproduces the metallicity dependence 
of the column density of the sharp transition from the atomic 
to fully molecular gas observed in the MW, LMC, and SMC. 

In order to investigate the environmental dependence of the 
star formation rate in the simulations, we perform a series of 
controlled test simulations. For each of these tests, we fix 
the dust-to-gas ratio in the H2 model and normalization of the 
emissivity of stellar particles at 1000 A to constant values and 
run the simulations for a significant period of time. 

We explore a grid of values of dust-to-gas ratio Dmw from 
10"^ to 1 .0 relative to the Milky Way value. The variable Z)mw 
scales the H2 on dust formation rate coefficient Rq and the 
absorption cross-section of dust in the Lyman- Werner band 
ctlw to the values characteristic for the Milky Way: 

o"lw = ^'mwO'o. (1) 
(Wolfire et al. 200 8) and 



Rd = Dmv/Ro 
3.5 X lO-i'' 



3 -1 
cm^ s 



where Rq 

o-Q = 2 X IQ -^^ cm^ ( |Draine & Be7toldi|1996nGlover & Mac 



LowpOOTa^ , respectively. 
The normalization of interstellar FUV flux at 1000 A: 



U^ 



= J 



ioooaZ-^mw 



used throughout this paper, is also defined to be in 
the units of the typical Milk y Way value Jmw - 



10" photons cm ^ s ' ster ' eV ' (|Draine|1978t|Mathis et al. 



[19831 . We explore the range of Umw frorn 0.1 to 100 in our 
test simulations. 

The star form ation model in our simulations closely follows 
the recipe 2 of Gnedin et al. ( 2009 | l with small numerical mod- 
ifications. Namely, the rate of star formation in each compu- 
tational cell with molecular fraction /h, > 0.1 is evaluated 
as 



dpi, 
dt 



TSF ' 



(2) 



where the time scale for star formation is defined as tsf - 
min(Tff , Tmax)- We follows the definition of [Krumholz & Tan| 
(|2007j) for the gas free-fall time. 



In 



(here p is the total mass density, including helium), and r^ax 
is the free-fall time in the gas with nsF = 50 cm"-*. We adopt 
esF - 0.005, which is lower than the value we adopted in 
Gnedin et al. ( 2009) and is still within the range of values 




Fig. 1 . — Average atomic-to-molecular gas transition as a function of total 
hydrogen number density for 9 test simulations (as distinguished by colors 
and line styles). 



advocated by Krumholz & Tan (2007 ). The lower value of esF 
that we adopt prov ides a better fit the THINGS measurements 
of the KS relation ( |Bigiel et al.|2008 1. 

The Tsf we adopt assumes that in low density cells, in 
which molecular fraction /h, is below unity, star formation 
proceeds mainly in unresolved molecular clouds on subgrid 
scales. This assumption then also motivates setting the maxi- 
mum free fall time to r^ax corresponding to the number den- 
sity of 50 cm"-* typical average density of molecular clouds. 
The /h, < 1 in these cells then can be viewed as reflecting the 
fraction of the total gas in such star forming molecular clouds, 
which themselves have /h^ = 1 , rather than incomplete con- 
version of the atomic gas into the molecular form inside the 
clouds. 

As we show below (see Fig.[7]and discussion in §[4]), the KS 
relation in our simulations is not very sensitive to variations of 
esF between 0.005 and 0.01 and nsF between 10 and 50 cm"^. 

3. THE ATOMIC-TO-MOLECULAR GAS TRANSITION 

The effect of two primary parameters, the dust-to-gas ratio 
Dmw ™d the interstellar FUV flux J/mw. on the transition 
from atomic to molecular gas is illustrated in Figure [T] as a 
function of the total hydrogen density, hh = "hi + "hii + 2nH, 
(the contribution of ionized gas «hii is negligible for densities 
shown in Figure [Th. As can be seen from the figure, both 
parameters affect the atomic-to-molecular transition in a non- 
trivial way. 

This scaling can be understood approximately if we ignore 
all physical processes except the formation of molecular hy- 
drogen on dust and dissociation of molecular hydrogen by the 
UV radiation in the Lyman- Werner band. This is necessarily 
an approximation, as many other processes are indeed impor- 
tant for the detailed balance of molecular hydrogen (see Ap- 
pendix), but the formation on dust and photo-dissociation are 
the dominant processes that control the atomic-to-molecular 
gas transition under normal ISM conditions. In this approxi- 
mation, the equilibrium abundance of molecular hydrogen can 
be determined from the balance of the formation and dissoci- 
ation rates (cf. Appendix) 



-o-uf/Nu 



RonHnm 



(3) 



where Flw = t^Mwro is the free space photo-destruction rate 
and Ro and ctlw are given by Equation (fill. The atomic gas 
becomes molecular only due to self-shielding and shielding 
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by dust (the last two factors on the left-hand-side of Equation 
([3]l). If the FUV flux is not too strong, the self-shielding by 
molecular hydrogen dominates; in this limit dust absorption 
can be neglected and Equation ^ becomes 

/h2 



-nn: 



(Equation (All i), so 



1 - /hj Umw ^oS h, 
where /hj = 1 "h and we ignore ionized gas. For our an satz 
for the self-shielding factor ^e, "h^^^ 
that 

All 

oc . 

1 - /h2 U mw 

Thus, the characteristic density at which molecular hydrogen 
fraction reaches a particular value (e.g., 50%) scales with the 
dust-to-gas ratio Dmw and the FUX radiation flux J/mw as 



MW 



MW 



4/7 



(4) 



In the opposite regime of large J/mw, the shielding by 
dust is expected to dominate over self-shielding, because self- 
shielding is a gradual function of the gas column density and 
may not be able to provide the required shielding for suffi- 
ciently large UV fluxes. In this regime. Equation Q becomes 



1-/h. 



-riii 



-e 



DmwO'o'Vh 



and the exponential factor is now large, so the characteristic 
column density for the atomic-to-molecular transition is 



ln(^Mw/^Mw) 



(5) 



Thus, as IGnedin et al. ( 2009 1 mention, in the regime where 
dust shielding dominates, the dependence of the characteristic 
column density on the FUV flux J/mw is only logarithmic. 

There is no way to convert between the characteristic col- 
umn density and the physical gas density easily. Neverthe- 
less, the following simple fitting formula captures the average 
dependence of the atomic-to-molecular transition on the dust- 
to-gas ratio and the FUV flux in our simulations: 



1 -I- exp {-Ax - 3x3) ' 



where x is given by 



x = A^/^ln(DMw^ 



(6) 



(7) 



Here n» = 25 cm 3, A is 

A = ln(l+gD^^(f/Mw/15)4/'), (8) 

and g is a fudge factor to approximately account for the tran- 
sition between the two regimes: g ~ \ when self-shielding 
dominates and g oc Dj^'^ when dust shielding dominates. 
We adopt the following fitting formula for the quantity g: 

\ + as + 



where 



0.04 



1 + s 



a - 5- 



1 + (f/Mw/2) 
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Fig. 2. — Average total hydrogen number density of atomic-to-molecular 
gas transition (defined as /h, = 0.5) as a function of the scaled dust-to-gas 
ratio Dmw and the FUV flux Uuw for all our test simulations. The point 
(Dmw = 0.001, 1/mw = 100) is missing because the resolution of our 
simulations is insufficient to capture the atomic-to-molecular transition in 
such extreme conditions. Solid lines show fitting formula of Equation |9}. 



and 



= 1.5 X 10"^ xln(l -h(3?7mw)' '') 



describes the transition to the regime when formation of H2 
via the gas phase reactions dominates. 

Figure [2] shows the value of the total (molecular, atomic, 
and ionized - although the contribution of ionized gas in all 
equations in this section is completely negligible) hydrogen 
density at which molecular fraction reaches /h, = 0.5 {x = 0). 
Our fitting formulae give the following approximate expres- 
sion for this density: 



A 



D 



(9) 



MW 



This equation is a better approxi mation than the the simple 
step-function ansatz proposed in ( [Gnedin et al.|2009] l. Figure 
|2] demonstrates that Equation |9] indeed provides an accurate 
model for the dependence of riuifu-, = 0.5) on Dmw and Umw- 

Figure [3] shows that Equation ^ works well for /h, > 
0.1 for alTsimulated cases (4 values of Umw and 7 values 
of Dmw), but it becomes somewhat less accurate for lower 
molecular fractions. The accuracy in the low /h, regime can 
be irnnroved with a simple modification: replacing x in Equa- 
tion ^ with x/g'^^. This change provides a more accurate fit 
for the range 10"^ < /h, < 0.1, but is less accurate than the 
above approximation for /h, > 0. 1 . Given that for modeling 
star formation the range /h, > 0. 1 is most relevant, we use the 
unmodified form of our fit as the fiducial approximation. 

Neither form of this fit describes the equilibrium H2 abun- 
dance (/h, ~ 10"^ - 10"**) in the Warm Interstellar Medium. 
Such a small abundance is, of course, not relevant to star for- 
mation. 

4. THE KENNICUTT-SCHMIDT RELATION AND ITS DEPENDENCE 
ON THE DUST-TO-GAS RATIO AND THE FUV FLUX 

The physics of the transition from atomic to molecular 
phase, discussed in the previous section, controls which local 
regions within the interstellar medium of simulated galaxies 
have high-molecular fraction and, hence, become the sites of 
star formation. Although the local rate of star formation in 
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Fig. 3. — Average atomic-to-molecular gas transition as a function of the fac- 
torized variable x (Equation JtJ) for all our test simulations (as distinguished 
by colors and line styles). The top panel shows the linear scaling of the y axis 
(most relevant for modeling star formation) while the bottom panel shows the 
y axis in log. Black squares on the right panel show the approximation from 
Equation |6|. 

these regions is sensitive to the parameters of the H2 forma- 
tion model and star formation recipe, the global star forma- 
tion rate surface density on larger, kiloparsec scales depends 
on the density and UV flux distribution within larger scales 
that are modeled self-consistently in the simulations. There- 
fore, once we fix the parameters of the model controlling the 
chemistry and star formation on small scales, we can exam- 
ine the predicted KS relation between the surface densities of 
various gas phases and the surface density of star formation 
averaged on large scale. 

Observationally, only the surface densities of atomic and 
molecular gas are directly measured and included in the es- 
timate of the "total" surface gas density. Eh- However, as 
we demonstrate below, the ionized gas may contribute signif- 
icantly to the total gas surface density under some conditions. 
Therefore, we deliberately avoid using the ambiguous nota- 
tion Sh and instead use the following notation explicitly indi- 
cating the components that are included in the surface density: 



-HI+HII+H; 



s E 



HII 



Ehi + E| 



«2, 



for the total surface density, uncluding both neutral and ion- 
ized gas, and 

for the surface density, including only neutral atomic and 
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Fig. 4. — Relation between Zsfr and the total surface density of gas (atomic, 
molecular, and ionized) for 9 different representative combinations of dust-to- 
gas ratio and the interstellar FUV flux (colored lines). The long-dashed line 
is the best fit relation of Kennicutt ( 1 998^ for z~0 galaxies. The gray shaded 
area shows the KS relation for the lo cal dwarf and normal spiral galaxies 
measured by the THINGS project ^Bigiel et al.|2008t . 



molecular gas. Note that we follow the observational prac- 
tice and do not include contribution of helium in the above 
gas surface densities. We emphasize again that in observa- 
tional work the total gas density is commonly identified with 
this second quantity. Eh - £hi+H2- 

As we mentioned in the previous section, this distinction 
is unnecessary for studying the atomic -to-molecular gas tran- 
sition on small scales, because the fraction of ionized gas is 
always small at densities at which the molecular fraction is 
significant. In other words, high-Znj regions are always sur- 
rounded by neutral atomic envelopes containing little ionized 
gas. However, regions of a kiloparsec scale can contain a 
mix of different ISM phases: from low-density ionized gas to 
high-density, molecular regions. In fact, diff'use ionized ISM 
gas is u biquitous in nearby galaxies (e.g., Hoopes & Walter- 
bos "2003). The warm (~ 10** K) diffuse ionized gas is present 
both inside the disk and at large distances (up to ~ 2 - 4 kpc) 
from the midplane both in the Milky Way (Reynolds 1989j 
I199H IGaensler et al. |l2008) and other nearb y galaxies (e.g., 
Hoop eseFal [|1999[ jC oflins et al. 2000; Ro ssa & Dettmar| 
2003 see Hattner etat.|2009| for review). This ionized gas can 
be a significant fraction of the total gas density. In the Milky 
Way, for example, the warm ionized gas accounts for ~ 25% 



of the total hydrogen col umn density of the disk (Reynolds 
1991[ Haffiier et al.|2009 l. One has to keep in mind the pos- 



sible presence of such gas in theoretical interpretations of the 
KS relation. 

For comparison with observations, the star formation rate in 
the simulations is averaged over 20 Myr and the gas and SFR 
surface densities are averaged on the scale of 500 pc. This 
specific choice corresponds to the averaging spatial scale and 
star formation indicator used in the T HINGS measurements 
( |Salim et al[][2007| [Bigiel et al.1[2008] l. We tested the sensi- 
tivity of the predicted KS relation to the specific choice of 
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the averaging temporal and spatial scales; such a comparison 
is presented in the Appendix (see Figure [T4|. Overall, the 
KS relation is robust to changes of spatial and temporal aver- 
aging scales with the range 0.5 - 2.0kpc and 20 - 100 Myr, 
respectively. Some modest trends are observed, but these are 
in general agreement with observations. 

In Figure |4] we show the relation between Ssfr and the to- 
tal surface density of gas (atomic, molecular, and ionized), 
2Hi+Hn+H2i for nine different representative combinations of 
dust-to-gas ratio and the interstellar FUV flux Dmw and C/mw- 
As could be expected, both the dust-to-gas ratio Dmw and the 
UV flux C/mw affect the relation significantly by affecting the 
atomic-to-molecular transition and the fraction of neutral gas 
in the ISM patches. Notably, the predicted Zsfr - 2hi+hii+h, 
relation does not agree with observations for any combination 
of Uym and Dmw- 

However, as we emphasized above, observational measure- 
ments often do not account for the contribution of ionized gas 
to surface density. We therefore present a separate predic- 
tion for the KS relation for the neutral gas only in Figure |5] 
for a representative subset of our test simulations. This figure 
demonstrates that the predicted Zsfr - ^^hi+h, relation for the 
parameter values representative of local galaxies (Dmw ~ 1 
and any value of Um w) is in good agre ement with both the 
older measurement of |Kennicutt| ( |1998) l and with the recent 
measurements by T he H I Nearby Galaxy Survey (THINGS) 
(Bigiel et al.|2008]l. In particular, our model approximately 



reproduces the rapid decrease of the SFR and increase of 
the scatter at Shi+Ht < lOMopc"^ and the change in the 
slope of the star formation rate vs gas surface density from 

SsFR ^Hi+H, ^° ^SFR 0^ ^et+H, ^HI+H, ~ 10^ MopC"^. 

The KS relations shown in Fig. Hand J5] can be accurately 
described by a simple fitting formula. Smce stars only form 
in molecular gas, the star formation rate surface density is 
proportional to the surface density of molecular gas, 



-SFR 



where tsf is the time scale for star formation (that may itself 
depend on the molecular gas surface density). If the neutral 
gas surface density Shi+h, is used as an argument, the reduced 
star formation rate at low gas surface density needs to be taken 
into account. 



1 



-SFR 



tsf(1 +S,/2hi+h,) 



(10) 



where is the characteristic surface density of neutral gas 
at which the relation steepens. At large gas surface densities 
(i.e., S» <& SH1+H2) we have: 



SsFR ~ (2^HI+Ho - ^m) ' 

Tsf 



(11) 



where is the saturation value of HI surface density, i.e. 
the maximum Ehi reached by gas as its total surface density 
increases to large values. Note that comparison of this equa- 
tion with the formula of Equation 10 shows that 



2 ■ 



Figure|5]demonstrates that, while the dust-to-gas ratio Dmw 
plays the dominant role in controlling the turnover in the 
£sFR-2^Hi+H2 relation at low surface densities for Dmw ^0.1, 



this is no longer the case at lower dust-to-gas ratios. Figure 
l6] shows the dependence of the characteristic "threshold" sur- 
face density S» on J/mw and Dmw for the full suite of our 
models. At Dmw < 0.1, 2* changes by an order of magnitude 
for J/mw changing by three orders of magnitude between 0.1 
and 100. Thus, although dependence of the KS relation on 
the FUV flux for higher dust content systems is expected to 
be weak, it can be stronger for dwarf galaxies at z a: and in 
high-z galaxies with low dust-to-gas ratios. 

The dependence of the H I saturation surface density on our 
two main parameters can be understood qualitatively if we as- 
sume that the density distribution in the ISM is approximately 
self-similar. Let us consider a large-scale region over which 
we measure the total hydrogen surface density Shi+hii+H2- 
Within this region the total hydrogen density has some density 
probability function (defined as a fraction of surface density 
contributed by gas of a given density «h), which in general 
depends on Zhi+hii+Ht: 0(«h> ^hi+hii+hJ- If the density dis- 
tribution is self-similar a region with a higher surface density 
will have more dense gas, i.e. 

0(nH,£HI+HII+H2) = >P(^), 

where ^ = nH/£Hi+Hii+H2 and 



f 

Jo 



The atomic hydrogen surface density is then simply 

2^Hi = f fm4>dn = £hi+hii+h, ( fm^(^)d^- 
Jo Jo 

If we assume that most of the atomic hydrogen mass is 
at densities near the atomic-to-molecular transition density 
«Hi^H2 (which is the case in our simulations), then we can 
use our parametrization from Equation (j6]) a function of fac- 
torized variable x (Equation (j7]i), so that 



nil 



dx 



dx Ehi+hii+h, A^V7' 



and 



fmn\i4'i^)dx. 



The last integral cannot be taken exactly, but given that the 
atomic -to-molecular transition is a rather steep function of the 
gas density, the integral can be approximated as 



1 



(/HinHlA)lHI^H, 



A3/7 
1 1 

■ nil i^H2 iA(^H i^H2 /2) Ax, 



(12) 



A3/7 2 

where Ax ~ 1 is the width of the atomic-to-molecular tran- 
sition (/hi - /h2 = 0.5) in the variable x, which should be 
essentially independent of any physical parameter. 

The saturation HI surface density ZJJj is obtained from 
Equation ( 12 1 in the limit of 2hi+hii+H2 — > oo, in which case 
the argument of if/ in Equation ( 12 1 can be replaced with zero, 
and we finally obtain 

A4/7 



»HI^H2 

2A3/7 



Ax oc 



D 



(13) 



MW 



We find that this scaling works well in our simulations, except 
in the Hmit of large Dmw and large J/mw, when the density of 
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Fig. 5. — KS relations for the neutral gas (atomic and molecular) predicted in models with different representative values for the dust-to-gas ratio (left panel) 
and the interstellar FUV radiation fluxe (right panel) are shown as colored lines. Dotted, short-dashed, and solid lines show the relation between Esfr and Shi, 
Shj, and Zhi+h, individually. The observed relations (long-dashed line and gray band) are the same as in Fig.|4] 




Fig. 6. — Characteristic threshold surface density Z, as a function of two 
main parameters Dmw and C/mw for all our test simulations. Cases with 
Dmw < 0.01 are not shown, as in our simulations gas at such low values of 
the dust-to-gas ratio never becomes fully molecular on 500 pc scale (and, 
thus, ZsF R can not be determined). The solid lines show the fitting formula of 
Equation {T4j. 




10 10= 



10= 



the ionized-to-atomic transition is not negligible compared to 
the density of the atomic-to-molecular transition. As a con- 
sequence, the contribution of the ionized gas is not negligible 
compared to the atomic gas, which leads to a decrease of 
compared to the value predicted by Equation (fT3]l. In the ex- 
treme case we consider (Dmw - 1, Umw — lOOjthe saturation 
H II surface density is 3-4 times higher than the saturation H I 
surface density. 

The following simple fitting formula corrects for this defi- 
ciency and provides a good fit for the characteristic "thresh- 
old" surface density, E«, and HI saturation surface density 



Fig. 7. — Dependence of the KS relation for the neutral gas (atomic and 
molecular) on the parameters of t he star formation recipe (2}. The long- 
dashed line is the best fit relation of lKennicuttH 1998) for z ~ Ogalaxies. The 
gray shaded area shows the KS relation tor the local dwarf and normal spiral 
galaxies measured by the THINGS project jBigiel et al.|2008) . 



= 2Et in all test cases we consider, 

1 



-2 



-y 1 + t^MWi^MW 



(14) 



The accuracy of this fitting formula is demonstrated in Fig- 
ures|6]and|9] For very low values of Dmw 5 0.01 the fit is not 
very accurate. This is most likely due to the limited volume 
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iDm«=0-3, U„^=10 
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Shi.h, (Me/pc^) 

Fig. 8. — Comparison of the KS relation for the neutral gas (atomic and 
molecular) for the full simulations and test runs which used Equation (6) to 
estimate the molecular fraction in the gas for a representative subset of values 
for £>MW ^nd J/mw- 

of our simulations: at such low dust-to-gas ratios the atomic- 
to-molecular transition shifts to extremely high gas densities, 
hh ~ 10"' cm""*, and our simulations lack 500 pc sized regions 
that would be dominated by such dense gas. Large volume 
simulations containing substantially more massive galaxies 
will be need to test the accuracy of the fitting formula ( 14 1 
in this regime. 

Finally, we have checked that our results are not particularly 
sensitive to the specific choice of the fiducial parameters esp 
and «SF- While the fiducial values provide the best fit to the 
median values of THINGS measurements (Bi giel et al. 2008) , 
a substantial variation in the adopted values for these parani- 
eters has only mild effect on our results, as we demonstrate in 
Figure [7] 

5. STAR FORMATION RECIPES 
5.1. Recipe for galaxy formation simulations 

In § [3] we have shown that atomic to molecular transition 
density can be well fit by fitting functions as a function of 
dust-to-gas ratio and FUV flux (e.g.. Equation (j6|). These 
fitting functions are an approximation to the average depen- 
dence of the molecular fraction on the total hydrogen density. 
The scatter in this relation around the mean may be important 
for particular observational measurements of the molecular 
abundance in the ISM. However, it is interesting to ask the 
question of whether we can reproduce results of our full sim- 
ulations by using the fit for molecular fraction given by Equa- 
tion (j6]l in star formation recipe of Equation (j2]i, instead of the 
true fii^ calculated using our full chemistry model. The results 
of such tests are shown in Figure [8] which demonstrates that 
using the fit to fn^inii) gives results closely matching results 
of the full calculations. 

This means that the approximation of Equation (j6]l can 
be used to implement the H2-based star formation recipe 
in galaxy formation simulations that do not follow the full 
molecular chemistry, provided that the resolution of the sim- 



ulations is sufficiently high (~ 100 pc) and that the values for 
the parameters Dmw and Umw could be estimated or assumed. 
The dust-to-gas ratio, Dmw can be estimated using local gas 
metallicity Z. Although the observed relation between Dmw 
and Z has a substantial scatter, on average the dust-to-gas ra- 
tio appears to be directly proportional to the gas metallicity. 



Zo' 



both for normal galaxies (Inoue 2003; Draine et al. 2007 



iOOTl 



Calura et al. 2008 ) and in low metallicity dw arfs OLisen: 
F'Ferrara 1998; Hirashita 1999l |Calura"eral . 2008; Madden 
2008||. Such a simple relation is, necessarily, a crude approxi 



mation, since not only the abundance, but even the properties 
of dust are known to be different in different galaxies. 

Relating the local FUV flux Umw is trickier, but sensible 
estimates can be made using the local SFR rat e averaged on 
a certain scale, as was done for example by IRobertson & 
Kravtsov] p008 j) . Given the steepness of the atomic to molec 



ular transition, the H2 -based star formation recipe amounts to 
the metallicity and FUV flux dependent density threshold for 
star formation. 

5.2. Star formation recipe for semi-analytic models 

The dependence of the KS relation on the dust-to-gas ra- 
tio and the FUV flux in our test simulations described in ^ 
can also be encapsulated by a simple recipe. Such a recipe 
can be used in semi-analytic models, in which radial depen- 
dence of gas surface density, star for mation, and chemical en- 
richme nt are modeled explicitly (e.g . , IFirma ni & Avila-Reese 
2000; Kravtsov et al.|2004||Dmtoire t al. 2007|^ 

As we discussed above, the dependence of the KS relation 
on Dmw and Umw in our models is due to the dependence of 
the characteristic HI surface density, E*, on these variables. 
We therefore parameterize the KS relation by the following 
fitting formula. 



-SFR 



2sfr(2hi+h2) 
(1+5:./Ehi+hJ' 



(15) 



where E» is given by Equation ( 14 1 and Esfr(2^hi+H2) is the 
star formation rate in the fully molecular gas at this surface 
density. For the latter, one can adopt either the original Ken- 
nicutt fit ( |Kennicutt|1998| ): 

y \1.4 



-SFR.K 



= 2.4X 10" 



kpc^ yr \ 1 Mq pc 2 



or the fit suggested by the study of |Bigiel et al 



2008] l: 



-SFR.B 



800 Myr 



max 1 



(16) 



(17) 



with the values of Sq, » lOOM^pc'^ and a x 0.5. Note that 
neither the slope at high surface densities a nor the character- 
istic surface density I.^ at which the slope steepens are well 
constrained by the current observations. 

Figure [9|sho ws that the fitting formula for the KS relation of 
Equation (|15[) together with Equation ([14]) reproduce the de- 
pendence oTthe KS relation on Umw aMZ^MW in simulations 
remarkably well. In semi-analytic models this formula can 
be used if one has some prescription for estimating Dmw and 
Umw in model galaxies. As we noted in the previous sections, 
these variables can be estimated approximately from the local 
metallicity of the gas and local star formation rate. 
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FiG. 9. — Scaled KS relation as a function of the neutral gas surface density, 
scaled by the characteristic surface density Z,. Black squares show the fitting 
formula jl5| with Z, given by Equation |l4f and Zsfr measured directly 
from the simulation as the star formation rate density in the molecular gas. 
Cases with Dmw < 0.01 are not shown, as in our simulations gas at such low 
values of the dust-to-gas ratio never becomes fully molecular on 500 pc scale 
(and, thus, £sFR cannot be measured). 



6. DISCUSSION AND CONCLUSIONS 

We have presented results of a phenomenological model for 
formation of molecular hydrogen and have illustrated the de- 
pendence of molecular fraction on the gas density, dust-to- 
gas ratio, and far UV radiation flux. We have also presented 
the large-scale Kennicutt-Schmidt relation arising in our sim- 
ulated galaxies when the local star formation is based on the 
density of molecular (rather than total) gas. Such approach al- 
lows us to avoid arbitrary density and temperature thresholds 
typically used in star formation recipes. Our results show that 
both the molecular fraction and the KS relation are sensititive 
to the dust-to-gas ratio and the FUV flux, although the sen- 
sitivity of the KS relation to the dust-to-gas ratio is stronger 
than to the FUV flux. 

We parameterize the dependencies observed in our simu- 
lations by fitting formulae (§ [3] and |4|i, which can be used 
to approximately account for H2 formation and H2-based star 
formation in simulations, which do not include a full H2 for- 
mation model and radiative transfer (see § |5.1| l. We demon- 
strate that our fitting formulae, when applied to realistic sim- 
ulations, produce results that are close to those obtained in 
simulations with the full H2 formation model and radiative 
transfer (Figure [Sll. 

We also provide fitting formulae for the dust-to-gas and 
the FUV radiation flux dependence of the KS relation that 
ca n be used in the semi-analytic models of galaxy formation 
(§ 5.2 1. One recent example of a model where such depen- 
dendcies can be relevant is the study of Dutton et al. ( 2009 ) 



density of atomic hydrogen is Shi ~ 10 Moyr"' and does 
not evolve with redshift. Our results, however, indicate that 
Shi should increase with increasing redshift, as metallicities 
(and, hence, the dust abundance) of galaxies decrease and 
their FUV fluxes increase. Conversely, the M»-MH,and SFR- 
M» relations should evolve differently if their expected depen- 
dence on the dust-to-gas ratio and the FUV flux is taken into 
account. Given that at lower metallicities (and, hence, the dust 
abundance) we expect smaller star formation rate for the same 
amount and spatial distribution of neutral gas, the trends de- 
scribed in this paper may potentially explain why the model 
of Dutton et al. (2009) overpredicts the specific star formation 
rate ( SSFR = SFR/M.) of small-mass galaxies at z > 3. 

One of the most interesting results of our simulations is that 
significant amounts of ionized gas can be present around high 
redshift gaseous disks. This ionized gas is akin to the d iffuse 
ionized gas observed in local galaxies (e.g., Hoopes & Walter- 
bbs 2003 i Haffner et al. 2009) and the Milky Way (Reynolds 
T589 1991 ;Gaensleretal. 2008). Our results indicate fliat the 
ionized gas may dominate the gas mass at low surface densi- 
ties (E < 10 Mq yr"'). Furthermore, our simulations show that 
ionized gas can remain a significant mass component at higher 
gas surface densities in environments with low dust content 
and/or high FUV fluxes (e.g., compare gas surface densities 
for a given Esfr in Figures|4]and|5]). One has to keep in mind 
the possible presence of significant amounts of ionized gas 
in theoretical interpretations of the KS relation and observa- 
tional estimates of the total gas mass. The significantly dif- 
ferent KS relation in the low dust-to-gas ratio, high FUV flux 
environments of high-redshift galaxies may also strongly bias 
gas mass estimates that use z = calibration of that relation 



(e.g., |Erb et al.|2006| |Mannucci et al.'2009|) 

As we discussed in Gnedin & Ki'avtsov| ( |2010| l, the dust- 
to-gas ratio and the FUV flux dependence of the KS relation 
that we observe in our simulations has a number of important 
implications for galaxy evolution, such as a lower efficiency 
of star formation in DLA systems, star formation confined to 
the highest gas surface densities of high-z disks, and generally 
longer gas consumption time scales in gaseous disks of high- 
redshift galaxies. The latter can be, at least partly, responsible 
for the prevalence of disk-dominated galaxies at low redshifts. 
This is because low efficiency of star formation can maintain 
disks gas rich until major mergers become rare. The outer, 
mostly gaseous regions of high-redshift disks should be more 
resistant against dynamical heating in mergers (e .g., |Robert- 
|son et al.|2004l [20061 [Springel & Hernquist|2005| ) and would 



help maintain f orming stella r disks dynamically cold during 
minor mergers (Moster et al. 2009) at later epochs. Moreover, 
minor mergers of forming disks should be largely gaseous, 
and gas brought in by such mergers should be deposited at 
large radii as it is ram pressure stripped by interaction with 
the gaseous disk and/or halo around it. This should prevent 
formation of large bulges, which was plaguing galaxy forma- 
tion models, and instead lead to formation of more extended, 
higher-angular momentum disks. This scenario is borne out 



in recent galax y formation simulations of Agertz, Teyssier, & 
Moore ( 2010[ l, who show that low efficiency of star forma- 



The results of that study indicate that the redshift evolution 
of SFR-M» relation of galaxies depends on the evolu tion of 
the re lation between stellar and molecular masses. DuttonJ 

|et al.| ( |2009. ) find that, in their model, the effective surface satellites can only be explained by a KS relation (Equation 



tion at high redshifts leads to more realistic disks and smaller 
bulge-to-disk ratios. 

Another interesting consequence of the complex depen- 
dence of the KS relation on the dust-to-gas ratio and th e FUV 
flux may be relevant to our own backyard. Recently, ( Orban 
et al.l200^ i noted that star formation histories of Milky Way 
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( 16 1) with the sharp threshold if the threshold varies semi- 
randomly within a modest dispersion of about 0. 1 dex. This 
variation is consistent with the variation given by Equation 
( [T4| for the values of Dmw and Umw typical for dwarf galax- 
ies (£>MW ^ 0.1, Umw ^ !)■ Since star formation histories 
of galactic satel lites are known to be highly variable (Ma^ 
|teo|[l998; Dolp hin et al.f2 005 ), the FUV flux is expected to 
vary accordingly; such variations may be responsible for the 
needed variation of the threshold in the KS relation, or, more 
precisely, the characteristic surface density from Equation 
(fT4li. 



20061 I Boissier et al.|2008||Wyder et al .120091 [Roychowdhury 



The high mass-to-light ratios (and hence low star formation 
efliciencies) of the Local Group dwarf spheroidal galaxies 
may also be partially explained by the environmental depen- 
dence of H2 abundance and, hence, star formation. Star for- 
mation in such low metallicity, low dust content dwarf galax- 
ies should be confined only to the highest gas surface densities 
(i.e., the central regions) while leaving the bulk of the gas at 
lower gas surface densities inert to star formation. This is con- 
sistent with observations of local dwarf low surface brightness 
galaxies which exhibit very low molecular gas fractions and 
anemic star formation rates ( [Matthews et al.|2005t [Das et aL] 



et al.|2009p . 

The examples described above illustrate the importance of 
further investigation of the effects of environmental depen- 
dencies of the KS relation discussed in this paper. The results 
and fitting formulae that we present should aid in implement- 
ing such dependencies in both cosmological simulations and 
semi-analytic models and should thus help to explore a wide 
range of possible effects. 
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APPENDIX 
H2 FORMATION MODEL 

In this Appendix we present the chemical reaction network of hydrogen and helium, as well as our phenomenological model 
for the formation of molecular hydrogen, in full detail (see also Gnedin et al. 2009 although we note that the model described 
here contains some modifications compared to the model used in this previous paper). 

We follow in detail 8 species of hydrogen and helium: H I, H II, He I, He II, He III, H2, H , and Hj . It is not, however, necessary 
to follow electrons separately, since, in all physical regimes of interest, abundances of Hj and H" are extremely small, so 

He ~ HhII + "He II + 2«HeIII- 

Note that this equation does not include any negative terms and thus «e will always be calculated with the relative error similar 
to the relative errors of neii, "He 11, and nneiii, but not larger. 

We follow all other species self-consistently and separately by solving the corresponding ODEs to avoid potentially unbounded 
increase of relative error in subtracting abundance of one specie from another (sometimes called "loss of precision"). For example, 
if the abundance of He III would be calculated by subtracting the abundance of He I and He II from the constant total abundance 
of He, the relative error of He III can be arbitratily large when the fraction of He III is small. 

We explicitly assume that all species are advected with the same peculiar gas velocity v. In this case the equations for the 
evolution of their number densities can be concisely represented as 



dn. 



1 



^ .• 3Hn: + -div,(n,v) = J,- + M: + £),•, 
at a ' 



(Al) 



where y = H I, H II, He I, He II, He III, H2, H , and Hj , the divergence is taken in comoving space x and three terms on the right 
hand side include reactions due to ionization balance, molecular chemistry, and dust chemistry respectively. This subdivision of 
the reactions into three sets is primarily for the sake of convenience and because we use different sources for different reaction 
rates. This separation is, of course, artificial - all the reactions take place together in a fluid element. 

The OTVET radiative transfer solver produces the radiation field at each computational cell that is used to calculate the rates 
for reactions between chemical species and radiation (including photo-ionization). We generically label these rates as F"^^ with 
various indicies. Since the self-shielding of molecular hydrogen and shielding by dust are not included in the OTVET solver, 
but are the ingredients of our empirical model, they are encapsulated into two factors, S h, and S d, with which we multiply the 
appropriate rates. These factors are described below. 

Ionization Balance 

Ionization balance terms include standard processes of photo-ionization, collisional ionization, and radiative recombination. 



and therefore only involve j - HI, HII, Hel, Hell, He III. We label all terms that include at least one of H2, H 
"molecular chemistry", and describe them all in the following subsection. 

' i'Hi = -"HirHi - Cmnenm + Rminenmi, 
i'Hii = -i'ei = -^Hii«enHii + "HirHi + Cmnenm, 

I Hel = -"HeirHel ~ CHeineWHel + (^^Hell + ^He Il)«e«He II> 

I Hell = -nHelirHell " (^^Hell + ^He Il)«e«He II " CHellWenHell + "HeirHel + C"Heine«HeI + ^He IIlW^nHe III> 
-^Helll = -■'^HemlelHelll + "HelirHell + CeelincnHelb 

Jh, = Xh- = Jh^ = 0. 



and HI as 



(A2) 
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Here Cj are collisional ionization rates, Rj are radiative recombination rates, a nd Dj are dielectronic recombination rates. For 
these rates we use highly accurate fitting formulae from Hui & Gnedin ( 1997 1. The recombination coefficients are computed 
self-consistently as a combination of case A and case B recombination, de]3ending on the gas opacity. 
The photo-ionization rates are derived from those returned by the radiative transfer solver and include the shielding by dust as 



Thi = SdTI\ [HI + 7^HII], 

Tnei = ^dF^i [Hel + r^HeH], (A3) 
rHeii = ^dF^u [Hell + He III]. 



In particular, we use the same factor to account for dust shielding in all three photo-ionization rates. Obviously, this is not exact, 
as the dust cross-section is a function of wavelength. However, since the effect of helium on molecular chemistry inside molecular 
clouds is thought to be small, helium ionization inside molecular clouds is sufficient to be treated rather approximately. 



Molecular Chemistry 

Molecular chemistry terms include a large set of reactions between H2, Hj , and H" and atomic species. The full set of equations 
we call "the full 8-species Model": 



Mm = FaHh- + Fb/ibj + 2FEnH2 + 2FLwnH2 - hnetim - hn-n-n-m - ^snHiiWHi - hn-n+n-ai- 

A;26nHeiinHi - 2^3onHi ~ '^h\n^-^n^, - 2^32nHi"Hei + 2ksnm\nn- + 2^6nenm + 'k-inn-,nmi+ 
Ikiiieriii^ + Ihtjiimnn. + 2^ionH2«H, + 2^iinHeinH2 + kunetiyi- + ^isnniWH- "+ ^2inmnH- + 
3A:22nH-nH+ + feneWHj + ^24nHeiinH2 + knnueinmi + ^28nHeiinH- + kignueinn-, 

Mmi = FBnH+ + 2FcnH: - ^snHiWHii - ^snH-WHii - kinu.nmi - ^leWH-nHii - hinHeinmi + k4nH*^nm+ 
^24nHeiinH2 + kienuinHsii, 

MhsI = -knnaunHei - kignH-nmi + fc24nHeiinH2 + ^25nHeiinH2 + ^26nHeiinHi + kimneiinH-, 

AIheII = -kiAnH.jiHeii - kunHjUHeii - ^26nHinHeii - kisnH-nHsU + ^27nHiinHei + fc29nH-nHei, 

AlHelll = 0, 

= -FdWh, - FsnHj - FlwWHj - kinH2nmi - kmenn. - hnmnn2 - ku)nn,nii^ - ^iinHeinH2- 
ki^rien^^ - ^24«HeiinH2 - ^snHeiiWH, + k2nn-nm + k4nn^_nm + k2in^nn- + fc3onHi+ 

^3inHl"H2 +^32nHl"HeI, 

Mh* = -FeKhj - FcKhj + FdHh, - k^nmnu* - k(,nen^ - fc2inH-nHj - ^22nH-nm + k3nninmi+ 

kiriii^nmi + ^lenHiiWH- + ^25nH2«Heii, 
Mh- = -FaHh- - kznmnn- - kstiminu- - k^nenu- - ^isWhiWh- - ^lenHiiWH- - kun^^n^— 

-kiinmnu- - ^28nHeiinH- - ^29«HeI«H- + kinetini + k23nenii2. 



(A4) 



where 



(■Fa = SuTf- [H-+r^HI + e], 

Fb = Sorf [H+-Hy^HI-HHII], 

Fc = Sd^I^ [H+ -hy ^ 2HII-He], 

' Fd = 5d5h,F^t [H2+y^H++e], 

Fe = 5d5h2F|T [H2+r^2HI (/iv> 13.6eV)], 

. Flw - '^D'5h2Flw [H2-t-y^2HI (Lyman- Werner band)]. 



The rate coefficients ki -ko,2 are taken from|Glover & Abel](|2008|; we do not list here all these reactions for bre vity. Cross section s 
for photo-rates A-D are given by Shapiro & Kang ( 1987), while the cross section for the reaction E is given by Abel et al. ( 1997 1, 
for both ortho- and para-H2. T he radiative transfer in the Lyman-Werner bands F^^ is treated fully self-consistently with 20,000 
frequency bins, as described in |Ricotti et al. (2002 1. 

Analogously to the previous section, we use the same 5'h2 factor to account for H2 self-shielding for reactions D, E, and LW. 
This is a crude approximation, but a more accurate treatment would introduce additional parameters that cannot yet be calibrated 
with the existing limited observational measurements. 

Equations ( A4i can be substantially simplified if we note that in all physical regimes relevant to cosmology the abundances of 
Hj and H" are always extremely small, so that they can always be assumed to be in the kinetic equilibrium, A1h+ ~ Mn- ~ 
(T. Abel, private communication). With this assumption and neglecting reactions involving k2\ and A:22, because their rates are 
oc nu-nii+ where both nn- and nu+ are small, expressions for the equilibrium abundances of H^ and H" can be derived in a closed 
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form, resulting in the following "6-species model": 



Ta + k2nm + k^nnii + k^ne + ^isWhi + ^leWHii + h^K^w + h^nn^i 
rD«H2 + hnmnmi + hriyi.nmi + ^lenHiinH- + ^25nH,nHeii 
Tb + Tc + k4nm + k^tie 

Mm = TAne- + rBnH+ + 2rEnH2 + 2rLwnH2 - ^in^wei - ^iWH-nHi - ^snHuHHi - hn^nm- 

kienHeiinm - 2^30«hi ~ 2^3i«hi"h, - 2^32nHi"Hei + 2^5nHiinH- + ^k^n^nu* + kjnu,niiii+ 
Ik^rien^, + Ikgriuinii, + Ikioriu^nn, + 2^nnHeinH2 + ^14«<.«h- + ^i5nHinH- + 

kunctlH^ + ^24nHeII«H2 + ^27nHeinHII + fc28nHeIinH- + ^29«HeI«H-, 

VHhii = rBnH+ + 2rcnH: - ^snHiWHii - ^5nH-nHii - ^vWHoHhii - ^lenn-nHii - fc27nHeinHii + fc4nH+nHi+ 

^24nHeIinH2 + ^nHlWHell, 
AIhcI = -^27«HIinHeI " ^29«H-nHeI + ^24«HeIinH2 + %«HeIinH, + ^26nHeIinHI + ^28nHeII«H-, 
AIhcII = -^24nH2nHeII - ^25«H2nHeII - ^26nHinHeII - kimH-nHeU + ^27nHIinHeI + ^29nH-nHeI, 

AIhciii = 0, 

A1h2 = -roWH, - rEnH2 - TlwWH, - ^7«H2«HII - ^8«e«H2 " hn^iTi^^ - ku)nn^n^^ - ^linHeinH2- 

kisnenn. - ^24nHeunH2 - ^snHeiiWH, + ^2nH-«Hi + ^4«mnHi + ^3onHi+ 

^3inHl"H2 +^32nHl"HeI- 

Finally, under normal ISM conditions the ionization balance of hydrogen and helium is controlled by the radiative recombina- 
tion, photo-ionization and ionization by cosmic rays. In this limit we can ignore all gas-phase molecular chemistry reactions, 

Mj ^ 0. 

We dub thi s approximation th e "minim al model" . The minimal model is often (justifiably) used in studies of local ISM (c.f. 
Kmmho lz & McKee 2005 ; Pelupessy et al. 2006 ; Krumholz & Tan 2007 ; Glover & Mac Low 2007a b ), but is also occasionally 
applied to high-redshift or low-metallicity systems (Krumh olz et al.|2 009a; Pelupessy & Papa dopoulos|2009| . We find, however, 
that the minimal model produces results that are reasonably close to the full model for £>mw <: 0. 1 (for any FUV flux), but 
becomes progressively less accurate for lower dust-to-gas ratios, mis-predicting the atomic-to-molecular transition as a function 
of density by a factor of 2 for Dmw ~ 0.01. 

In order to maintain high accuracy for the full sampled range of Dmw ™d Um^n, oil simulations presented in this paper were 
performed with the 6-species model. 

Dust Chemistry 

In our model the only dust chemistry reaction that we include is the formation of molecular hydrogen on dust. 



•Oh, = DMw^()CpnHi(nHi + 2nHj, 
£)hi = -2i)H2, . 

2)hII - ^Hel - f Hell = ^^Helll = ^'h" 



(A7) 



where Rq - 3.5 X 10"'^ cm^ s"' ( [Wolfire et al.|2008 see Equation (jl]l) and Cp is the clumping factor inside molecular clouds 
which takes into account the fact that the gas is clumjjed on subgrid scales unresolved in our simulations (also see lGnedin et al. 
2009 1. The clumping factor Cp is a parameter of our model, we discuss a reasonable choice for its value below, in 5 A. 6 



Heating, cooling, and thermodynamics 

For the heating and cooling terms in the equation for the internal energy we include all of the terms normally included in the 
simulations of first stars and in the ISM models. Specifically, the entropy term in the energy equation for the gas can be written 

as 

dt 

where "/V and C are heating and cooling terms, 

C - Cci + Crr + Cder + Cle.a + Cff + Cqx + Cle,H2 + Cle,z + Co- (A8) 
In the heating function, we include 



"Hpi: : photoionization heating due to H I, He I, and He 11, using cross-sections from Hui & Gnedin ( 1997 1; 



'HcMB- ■ Compton heating/cooling on the CMB (Hui & Gnedin 1997 1; 



'HLya- ■ heating by Lya photons ( Tozzi et al. 2000 1; 
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Fig. 10. — Cooling functions (per hydrogen nuclesus) for 3 representative values of gas metallicity Z and the FUV flux Umw- In this plot we assume 
Dmw = Z/Zq. Blue points show the full cooling function (including all relevant physical processes), while red points show the result of excluding H2 cooling. 
Black lines trace the H2 cooling function from ^Galli & Pallaj. 1998.1 (left panel) and the standard, metal-free cooling function (right panel). 



•TYh,: : heating due to photo-dissociation of H2, TYh, = 0.4 eV x rin^ (Fd + Fe + Flw) (Equation (|A5|l); 
"J^PAH- ■ photo-electric heating on PAH, implemented as in ( 



Glover & Mac Low 



( 2007a 1; 



'HcR- '■ cosmic rate heating, assuming that the cosmic rate density scales as the dust-to-gas ratio, implemented as in Glover & 
I Mac Lowlp007a| ). 

CooUng processes include 



Cci- '■ cooling due to collisional ionizations of HI, He I, and He 11 (Hui & Gnedin 1997 1; 



Crr: : cooling due to radiative recombinations of HII, Hell, and He III (Hui & Gnedin 1997 1; 



Cder: : cooling due to di-electronic recombination of He III (Hui & Gnedin 1997 1; 



Cle,a'- ■ line exitation cooling of H I and He II ( Hui & Gnedin 1997 1; 



Cpp: : free-free emission ( Hui & Gnedin 1997 1; 

Cqx: : cooling due to charge exchange reactions between H2, H , HI and free electrons (reactions 8, 9, 10, 14, and 15 from 
[Glover & Abel (2008)); 



Cle,H2: : line exitation coohng of H2 (Glover & Abel 2008 1; 



Cle.z: : line exit ation cooling of heavy elements, using Sutherland & Dopita ( 1993 1 cooling functions for T > 10 K and Penston 
(|1970|) and|Dalgarno & McCrayl(|l972|) rates in the T < 10"* K regime; 



Cd- '■ cooling on dust from 



Draine 



(1981 



Some of the reaction rates involving H2 depend on the ortho-to-para ratio of molecular hydrogen. For this ratio and other 
thermodynamic quantities {j{T), U{T), etc) we use exact expressions computed from quantum-mechanical statistical sums (Turk 
et al, 2010, in preparation). 

Examples of cooling functions from our simulations are given in Figure 10 The cooling function, in general, is not a function 
of gas temperature only, but also depends on the gas metallicity Z, the energy density of the incident radiation field t/y, the 
number density of baryons ni, (although for rtt, < 10"^ cm"^ the dependence on the last two parameters always enters as Uv/rib), 
and abundances of all atomic and molecular species Xj = rij/rit,. Therefore, when plotted as a function of temperature, the cooling 
function takes a range of values (depending on the values of other gas properties) rather than a single, unique value. 
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Interestingly, Figure 10 shows that the cooling rate at T < lO'^K is dominated by cooling due to molecular hydroge n, rather 
than by low ionization metal speci es such as OI or CII. Molecular hydrogen coohng i s often assumed to be negligible (c.f. Wolfire 
et al.|2 003', 'Stabler & Palla'2005') due to lower cooling rates (c.f.lGalli & Palla(1998) 
rates of Glover & Abel (2008), which are considerably higher thi 



However, we use the updated H2 cooling 
an the previous estimates. As Figure [TO] shows, the new H2 



considerably 

cooling rates dominate over the low ionization metal species at T < 5000 K 

Shielding Factors 

The two shielding fac tors. So and Sh2, together with the clumping factor Cp, are important parameters of our empirical model. 



(A9) 



As Gnedin et al. ( 2009 | l explain, we use an ansatz similar in spirit to the Sobolev approximation to estimate dust shielding: 

- g-Duwo-oinm + 2«H2)i'Sob 

where Dmw is the dust-to-gas ratio in units of its Milky Way value (see §|2]i, ctq = 2 x 10"^^ cm^, and 

isob = p/(2|Vp|). (AlO) 

Note that the value for o-q that we use in this paper is twice lower than the one listed in Gnedin et al.'('2009V the new value is 
a commonly adopted value for this parameter for the Milky Way type dust, and provides a better quantitative fit t o the existing 
observ ational constraints. In addition, a factor of 2 in the denominator of the expression for Lsob was missing in [Gnedin et al. 
(|2009|) - this was a typo, and the correct expression was used when simul ations were ru n. 



TEe major change bet ween our curre nt model and the model of Gnedin et al. ( 2009[ l is in the form of the molecular hydrogen 
self-shielding factor. In Gnedin et al. (2009 ) this form was modified from the commonly used formula of Draine & Bertoldi 
( [T996) , because the FUV flu x in Gnedin et al. ( 2009 ) w as much higher than the Draine value. In our present tests, we find that 
we can use either the original Draine & Bertoldr| ( |1996| l formula or their simpler and more approximate expression. 



(A?H2/10''*cm-2) ^''^ , forA^H, > lO'^^cm^^ 



for A^H, < lO'^cm-^, 



(All) 



which we actually use for computational efficiencMj 

Finally, to complete the full specification of our ch emica l model, we need to estimate the column density of the molecular gas, 
Nu^, for the self-shielding factor given by Equation (All 1. Unfortunately, we cannot simply use the Sobolev approximation to 



derive A^h, similar to the column density of dust in Equation ( A9 1, because H2 absorption is concentrated in separate absorption 



lines and is sensitive to the internal velocity dispersion inside molecular clouds. These velocities are unresolved in our simula- 
tions, but can greatly reduce the self-shielding of molecular gas. Dust, on the other hand, absorbs UV radiation in continuum and 
is thus not affected by velocity distribution of the gas. 



Therefore, we introduce the following simple ansatz for the effective column density Na^ for Equation (All 



(A12) 

where Lc is the velocity coherence length of the molecular hydrogen inside molecular clouds. Since we cannot deduce this 
quantity from observations or other calculations, we treat it as another parameter of our model. 

With the expressions for the shielding factors above, the only two parameters of our model are Cp and Lc- These parameters 
can only be determined by comparing the simulation results to the observational data. 

Calibration 

As the primary data sets used to calibrate th e model, we use the measurements of atomic and molecular gas surface densities 
in nearby spirals from Wong & Bhtz ( 2002b) and measurements of gas fractions along the lines of sight to individual st ars for 
atomic (Goldsmith & Li 2005 ) and molecular gas in the Milky Way and Magellanic Clouds (Tumlinson et al.,2002 ; Gillmon et al. 
[2006; Wolfire et al. 2008). 

We calibrate the two parameters of the model: the clumping factor Cp and the molecular coherence length Lc. We find, however, 
that there is no unique best-fit set of parameters. Instead, any combination of these two parameters that satisfy the constraint 



LrCr, 



10 pc 



provides an acceptable fit to the observational constraints. As an example, we show on the left panel of Figure 1 1 fits to the I Wong 
[& Blitz_(2002b) measurements (averaged over all galaxies they observed) for three combinations of the parameters and Cp. In 
general, higher clumping factors result in the lower atomic contents at high surface densities, but the trend is too weak to be of 
any statistically significant constraining power. 

As a fiducial set of parameters we choose the combination Lc = 0.3 pc and Cp - 30. This choice provides a marginally better 
overall fit to the obse rvations, and is also consistent with estimates of the gas clumping factor deep inside molecular clouds 
( McKee & Ostriker|2007^ ). The fiducial value of Cp is somewhat larger than the estimates of the clumping factor from numerical 

simulations of turbulent molecular clouds, Cp - e^^^p, where o"inp ~ 1 - 1 . 5 is the dispers i on of the lognormal density distribution 



inside the clouds. However, the value of Cp - 10, which was used in IGnedin et al.l (12009') and is more consistent with the 



numerical simulations of turbulent molecular clouds would provide an almost equally good to the existing observations, if it is 
used with Lc ~ 1 pc. 



^ We have indeed verified that a more complex formula (Equation (37) of 
[Draine & Bertoldi 1 1996 1) produces essentially indistinguishable results from 



the more approximate form of Equation |A1 1 
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Fig. 11. — Average atomic and molecular gas surface densities as functions of tlie total (neutral) hydrogen gas surface density averaged over 500 pc scale. The 
left panel show three test simulations with three values of the clumping factor Cp and molecular coherence length L^. Filled squares and open circles with error 
bars mark the observed average atomic and molecular hydrogen surface densities at = 10, 30, and 100 A/© pc"- from Wong & Blitz 2002b| . The right panel 
shows our fiducial model (Lc = 0.3 pc, Cp = 30) together with the rms scatter (shaded bands) around the averages. The error-Dars on me oDservational points 
now show the dispersion around the average rather than the error of the mean. 
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Fig. 12. — Atomic (bottom) and molecular (top) gas fractions as functions of the total (neutral) hydrogen gas column density along individual lines of sight 
through the galactic disks. Colored points shows our fiducial test simulation (L^. = 0.3 pc, Cp = 30), while black points show observational measureme nts. The 
left panel shows the (Dmw = 1, t^MW = 1) simulation case and the observational measurements of molecular fractions in the Milky Way galaxy from Gillmon| 
[et al. 1 2006 1 (filled triangles) and Wolfire et al. ( 2008 1 (filled squares) and atomic fractions measurements from Goldsmith & Li 1 2005 1. The right panel shows 
(Dmw = 0.3, Umw = 10) (blue points) and (L*mw = 0.1, J/mw = 100) (red points) simulation cases that should bracket possible values of these para meters for 
Magellanic Clouds. Filled saquares and triangles on the top panel show the measurements for LMC and SMC molecular fractions respectively (Tumlinson et al.| 
|2002^ . On the bottom panel the measurements are for SMC ^Leroy et al.|2007j l, to be compared with red points. 



Dependence on Numerical Resolution 

Any sub-cell model would be of limited value, if it was only applicable to a narrow range of numerical resolutions. In order 
to test the range of spatial resolutions over which our model performs robustly, we have re-run a subset of our test simulations, 
varying the maximum allowed level of refinement between 6 and 10, compared to our fiducial value of 9 (cell size of Ax - 65 pc 
at z = 3 in physical units). 

The results of these tests are shown in Figure 13 for the atomic-to-molecular transition and the KS relation. In order to perform 
a genuine resolution test, in each run with different resolution we only show cells that are refined to the lowest allowed level. For 
example, in the run with the maximum level 10, we only show cells from level 10, so that level 9 cells, which are also present 
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Fig. 13. — Dependence of the atomic-to-molecular transition (left) and the KS relation (right) on numerical resolution in our model. The left panel shows three 
representativel cases (Dmw, Umw) = (U !)• (Dmw, Umw) = (01, 100), and (Dmw, Umw) = (0.01, 100), while only the first case (Milky Way like parameters) is 
shown on the right panel for the sake of clarity (the other two cases show similar behavior). The value of the cell size Ax on the highest resolved level is shown 
for each line. Black squares on the left panel trace the approximate fit j6}. 

in that test run, do not contaminate Fig. |T3] Of course, in realistic simulations cells from all levels that contain molecular gas 
are going to contribute to the /h, - «h relation, so Fig. 13 actually exaggerates the effect of changing resolution. At resolutions 
Ax < 260 pc our model performs robustly down to the smallest scales we are able to probe (Ax x 30 pc). At coarser resolution of 
Ax = 520 pc small molecular clouds in low density gas are not captured properly, resulting in a sharper fall- off in the KS relation 
at low values of Shi+Hj • In addition, the Sobolev-like approximation for the dust column density (Equation ( AlO i) overestimates 
the column density significantly, which results in the atomic-to-molecular transition shifting towards lower density gas (especially 
for low dust-to-gas ratio and high FUV flux). We conclude, therefore, that spatial resolution of at least 250 pc is required for our 
model to work robustly. 



Dependence on Averaging Scales 

The exact value of the star formation rate surface density and the gas surface density in principle can depend on the specific 
choices for the spatial and temporal scales over which Eh and Ssfr are averaged. Observational studies (Kennicutt 1998 ; SalimJ 
et al.|2007 |Bigiel et al.|2008) often use a combination of star formation estimators that correspond to different temporal scales. 



Therefore, the best approach would be to model the observational methodology exactly, but this is not feasible in practice. In this 
paper we adopt a simplified procedure, and select the fixed values for both the temporal Af and spatial AZ averaging scales. The 
sensitivity of our results to the exact choice for these two scales is shown in Figure[T4] In general, the KS relations measured in 
the simulations are robust for Af < 30 Myr and A/ < 1 kpc. For larger spatial and temporal scales modest trends are observed. 
Several processes can contribute to such trends. For example, if the star formation at low surface densities is intermittent on the 
time scale of the averaging (i.e. stars form only during episods of duration comparable to the averaging time period), the average 
SsFR can depend on the time peri od used for averagin g. This may explain the weak trend at low Eh with Af. Such trend is also 
consistent with observations (e.g., |Boissier et al.|2007| l, which show that star formation derived from the UV flux is more spatially 
extended compared to the star formation derived from H„, which corresponds to time period of ~ 10^ years. Overall, our results 
are quite robust to changes of spatial and temporal averaging scales within the range of values used in observations. This relative 
insensitivity of the KS relation (besides the weak trends mentioned above) is in general agreement with observations, which 
indicate bro adly consistent KS relations derived usi ng different star formation indicators and a wide range of spatial averaging 
scales (e.g., [Kennicutt et al.|2007||Bigiel et al.|2008| l. 
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